clear all
use "datasets/MRA_waste_sites"

*In all cases, the WLS-RE model with reg command without clustered errors is used. The unclustered error calculation is necessary for "cv_regress" command and not harmful for BT error calculation as only the coefficents are of importance and not the standard errors. In all cases aweights are used for the ease of coding instead of manually weighted observations. This does not alter the resulting coefficients.

*creates the x-axis for the scatter graphs to detect outliers
sort elas
gen obsno_elas = _n

*----------------------------------------------------------------------------*
* Out-of-sample error based on meta-regression results
*----------------------------------------------------------------------------*

quietly metareg elas elas_var publish year_publish_c site_m NPL_n NPL_y NPL_u active_n active_u job_n job_u non_haz nuclear air water element_u North_America Asia Other_continent cleanup0 cleanup2 cleanup3 cleanup_unclear dist_greater_mean sample_c sale_ind num_expl_c oth_disamen oth_amen access industry demoecon time_control direction interaction log_log OLS_spatial, wsse(elas_SE) tau2test

scalar tau2 = e(tau2)
display tau2
gen elas_var_re_pb_sq= elas_var + tau2
gen precision_sq_re_pb_sq= 1/elas_var_re_pb_sq

quietly reg elas elas_var publish year_publish_c site_m NPL_n NPL_y NPL_u active_n active_u job_n job_u non_haz nuclear air water element_u North_America Asia Other_continent cleanup0 cleanup2 cleanup3 cleanup_unclear dist_greater_mean sample_c sale_ind num_expl_c oth_disamen oth_amen access industry demoecon time_control direction interaction log_log OLS_spatial [aw=precision_sq_re_pb]

cv_regress, genhat(elas_p_out)
gen MRA_full= abs((elas-elas_p_out)/elas*100 )
sum MRA_full, detail
local max=r(max)
gen MRA_full_wo_max= MRA_full if MRA_full <`max'
sum MRA_full_wo_max, detail // reduces mean from 1602 to 543
/* omitting Walsh and Mui (2017) study that represents the vast majority of extreme APE
gen MRA_full_wo_max_88 = MRA_full_wo_max if ID_Study !=88
sum MRA_full_wo_max_88, detail // reduces from 543 to 213, dropping 14 obs
*/

*Non-hazardous
quietly metareg elas elas_SE publish year_publish_c site_m  active_n active_u job_n job_u  air water element_u North_America Asia Other_continent dist_greater_mean sample_c sale_ind num_expl_c oth_disamen oth_amen access industry demoecon time_control direction interaction log_log OLS_spatial if non_haz==1, wsse(elas_SE) tau2test
scalar tau2 = e(tau2)
display tau2
gen elas_var_non_haz= elas_var + tau2
gen precision_sq_non_haz= 1/elas_var_non_haz

quietly reg elas elas_SE publish year_publish_c site_m  active_n active_u job_n job_u  air water element_u North_America Asia Other_continent dist_greater_mean sample_c sale_ind num_expl_c oth_disamen oth_amen access industry demoecon time_control direction interaction log_log OLS_spatial if non_haz==1  [aw=precision_sq_non_haz]

cv_regress, genhat(elas_p_outs_non_haz)
gen MRA_non_haz= abs((elas-elas_p_outs_non_haz)/elas*100 )
*scatter MRA_non_haz obsno_elas,mlabel (ID_Study) // one observation
sum MRA_non_haz, detail
local max=r(max)
gen MRA_non_haz_wo_max= MRA_non_haz if MRA_non_haz <`max'
sum MRA_non_haz_wo_max, detail // reduces mean from 13426 to 463
*scatter MRA_non_haz_wo_max obsno_elas,mlabel (ID_Study) // 

*NPL
quietly metareg elas elas_var publish year_publish_c site_m NPL_n NPL_y NPL_u active_n active_u job_n job_u non_haz nuclear air water element_u North_America Asia Other_continent cleanup0 cleanup2 cleanup3 cleanup_unclear dist_greater_mean sample_c sale_ind num_expl_c oth_disamen oth_amen access industry demoecon time_control direction interaction log_log OLS_spatial if NPL==1, wsse(elas_SE) tau2test

scalar tau2 = e(tau2)
display tau2
gen elas_var_NPL= elas_var + tau2
gen precision_sq_NPL= 1/elas_var_NPL

quietly reg elas elas_var publish year_publish_c site_m NPL_n NPL_y NPL_u active_n active_u job_n job_u non_haz nuclear air water element_u North_America Asia Other_continent cleanup0 cleanup2 cleanup3 cleanup_unclear dist_greater_mean sample_c sale_ind num_expl_c oth_disamen oth_amen access industry demoecon time_control direction interaction log_log OLS_spatial if NPL==1 [aw=precision_sq_NPL]

cv_regress, genhat(elas_p_outs_NPL)
gen MRA_NPL= abs((elas-elas_p_outs_NPL)/elas*100 )
*scatter MRA_NPL obsno_elas,mlabel (ID_Uni) // 1 obs sticks out
sum MRA_NPL, detail
local max=r(max)

/* 
*gen MRA_NPL_wo_max = MRA_NPL if MRA_NPL <`max' 
*sum MRA_NPL_wo_max MRA_NPL, detail // reduces mean from 359 to 202

*omitting Walsh and Mui (2017) study that represents the vast majority of extreme APE
scatter MRA_NPL_wo_max obsno_elas,mlabel (ID_Uni) // 3 obs stick out
gen MRA_NPL_reduced = MRA_NPL if MRA_NPL <=2500
sum MRA_NPL_reduced , detail // reduces mean from 359 to 125
*/
*---------------------------------------------------------------*
*Out-of-sample errors based on Unit value transfer
*---------------------------------------------------------------*
*full sample
quietly metareg elas, wsse(elas_SE) tau2test
scalar tau2 = e(tau2)
display tau2
gen elas_var_re_full= elas_var + tau2
gen precision_VT_full = 1/(elas_var_re_full)

quietly reg elas [aw=precision_VT_full]
cv_regress, genhat(elas_p_out_VT)

gen VT_full= abs((elas-elas_p_out_VT)/elas*100 )
sum VT_full, detail
local max=r(max)
gen VT_full_wo_max= VT_full if VT_full <`max'
sum VT_full_wo_max, detail // reduces mean from 5247 to 681
/*
gen VT_full_wo_88 = VT_full_wo_max if ID_Study !=88
sum VT_full_wo_88 // reduces from 679 to 219, dropping 14 obs
scatter VT_full_wo_88 ID_Study, mlabel(ID_Uni)
*/

*Non-haz
quietly metareg elas if non_haz==1, wsse(elas_SE) tau2test
scalar tau2 = e(tau2)
display tau2
gen elas_var_re_non_haz= elas_var + tau2
gen precision_VT_non_haz = 1/(elas_var_re_non_haz)

quietly reg elas if non_haz==1 [aw=precision_VT_non_haz]
cv_regress, genhat(elas_p_out_VT_non_haz)
gen VT_non_haz= abs((elas-elas_p_out_VT_non_haz)/elas*100 )

*scatter VT_non_haz ID_Study, mlabel(ID_Study)
sum VT_non_haz, detail
local max=r(max)
gen VT_non_haz_wo_max = VT_non_haz if VT_non_haz < `max'
sum VT_non_haz_wo_max, detail // reduces from 39780 to 601, dropping 14 obs
*scatter VT_non_haz_wo_max ID_Study, mlabel(ID_Uni)

*NPL
quietly metareg elas if NPL==1, wsse(elas_SE) tau2test
scalar tau2 = e(tau2)
display tau2
gen elas_var_re_NPL= elas_var + tau2
gen precision_VT_NPL = 1/(elas_var_re_NPL)

quietly reg elas if NPL==1 [aw=precision_VT_NPL]
cv_regress, genhat(elas_p_out_VT_NPL)
gen VT_NPL= abs((elas-elas_p_out_VT_NPL)/elas*100 )

*scatter VT_NPL ID_Study, mlabel(ID_Uni)
sum VT_NPL, detail 
/*
gen VT_NPL_wo_max = VT_NPL if VT_NPL <`max'
sum VT_NPL_wo_max, detail // reduces from 133 to 
scatter VT_NPL_wo_max ID_Study, mlabel(ID_Uni)
*/

*--------------------------------------------------------
* Table 5: Benefit transfer errors
*--------------------------------------------------------

* with full sample 
quietly estpost tabstat MRA_full_wo_max MRA_non_haz_wo_max MRA_NPL VT_full_wo_max VT_non_haz_wo_max VT_NPL , column(statistics) statistics(n mean sd median min max)
ereturn list
esttab using "tables/Table 5. Benefit transfer.rtf", replace cells("count mean(fmt(2)) sd(fmt(2)) p50(fmt(2)) min(fmt(2)) max(fmt(0))")

*----------------------------------------------------------*
*Figure A7: Absolute percentage transfer error
*----------------------------------------------------------*

sort elas
twoway line MRA_full_wo_max obsno_elas, xtitle("Number of observations") ytitle("Absolute percentage transfer error")
graph rename Graph MRA_full, replace
graph export "figures/Figure A7. Absolute percentage transfer error.tif", replace as(tif)


clear


